Fitting Skyrme functionals using linear response theory 
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Abstract 

Recently, it has been recently shown that the linear response theory 
in symmetric nuclear matter can be used as a tool to detect finite 
size instabilities for different Skyrme functionals U[2]. In particular 
it has been shown that there is a correlation between the density 
at which instabilities occur in infinite matter and the instabilities in 
finite nuclei. 

In this article we present a new fitting protocol that uses this 
correlation to add new additional constraint in Symmetric Infinite 
Nuclear Matter in order to ensure the stability of finite nuclei against 
matter fiuctuation in all spin and isospin channels. As an applica- 
tion, we give the parameters set for a new Skyrme functional which 
includes central and spin-orbit parts and which is free from instabil- 
ities by construction. 



1. Introduction 

The formalism of the Linear Response (LR) theory has 
been presented in a recent series of articles written by 
the NESQ collaboration [31 HI [S] for the case of a general 
Skyrme functional |6] including spin-orbit and tensor in- 
teractions. Such a formalism can be adopted as a useful 
tool to detect instabilities that plague Skyrme function- 
als [H IIS]. 

In the present article, we present an improved fitting 
protocol to build Skyrme functionals, that, for the first 
time, includes the stability of Symmetric Nuclear Mat- 
ter (SNM) up to a given density as an additional con- 
straint. Such method is more powerful than the standard 
constraints on the Landau parameters, since the latter are 
only valid in the so called long-wavelength limit and are 
therefore not able to detect finite-size instabilities. 

An example of this kind of instabilities has been revealed 
by Lesinski et al. [8j and was shown to be related to large 
values of C^^ , the coupling constant of the term piApi 
in the Skyrme functional [H]- This coupling constant, if 
not properly constrained during the fitting procedure, can 
drive the merit function used to adjust the parameters 
of the interaction in a region of matter fluctuations that 
give unphysical properties for the ground state of some 
nuclei. Since this gradient term does not contribute to 
the Landau parameters, this kind of instabilities is only 
detectable trough the LR method. 
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The article is organized as follows: in Section [5] we re- 
call the main properties of a Skyrme functional and in Sec- 
tion[3]we briefly remind the LR formalism; in Section|4]we 
present the results of our new fitting method and finally 
we give our conclusions in Section [5] 

2. The energy functional 

The density-dependent Skyrme interaction considered in 
the present article has the standard form 

+ i<3(l + a;3P.)Po(R)'5(r) 



^h{l + xiP,) [k'2<5(r) + J(r)k2 
t2{l + X2pa)yi' ■ S{r)k 
-iWo(<Ti +0-2) • [k' X S{r)k] , 



(1) 



where we used r = ri — r2 and R — ^ (ri -|-r2) for 
the relative and center of mass coordinates, respectively. 
Wc also defined Pc the spin exchange operator, while 
k = — ^(Vi — V2) is the relative momentum acting on 
the right and k' its complex conjugate acting on the left. 
Finally po (R) is the scalar- isoscalar density of the system. 
The Energy Density Functional (EDF) is derived by keep- 
ing all the terms except the Coulomb exchange one which 
is treated using the Slater approximation. The center of 
mass correction is approximated by its one-body part. 

3. Linear response theory 

The Random Phase Approximation (RPA) response func- 
tion in each channel (a) — {S, M, I) = (spin, projec- 
tion of the spin and isospin) is obtained by solving the 
Bethe-Salpeter equations for the correlated Green func- 
tions G^^pj^ written as 

G'gpA(g,w,ki) = GHF('?,i:^,ki) + GHF(g,w,ki) 



(27r 



^3l7^^"(q,ki,k2)G^^PA(<Z,C.,k2), (2) 



where Ghf('Z, w, ki) is the unperturbed particle-hole prop- 
agator, and Vp^" (g,ki,k2) is the residual interaction ob- 
tained through second functional derivative techniques. 
We refer to articles [71E1I1] for more details. Once Eqs. (0) 
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Fig. 1: (Colors online). In the left panel, we show the position of the poles in SNM for the SKP functionals. The horizontal 
black dotted-dashed line represents the saturation density of the system, psat, while the red dotted line is placed at 1.4psat. In 
the right panel, we show the response function in SNM for the value of density and transferred momentum indicated by the 
black arrow in left panel. See the text for details. 



are solved, we can calculate the response function of the 
system with 

x"(^,9)=4y'^Gg,pA(9,c.,k), (3) 

where uj is the transferred energy and q is the transferred 
momentum. As usual we chose the momentum q oriented 
along the z— axis of the reference frame. We also use nat- 
ural units ?i = c = 1 to simplify the notations. When 
an instability appears in the infinite system, it manifests 
itself by a pole in the response function x"(cLi,g) at zero 
energy. 

In Figure [1] we show the positions of the poles in SNM 
for the Skyrme SkP [TU] functional. We observe that 
in the channel {S,M,I) — (0,0,1), there is a pole for 
Ppoig = 0.205 fm~^ and q — 2.88 fm^^. Lesinski et al. [5] 
showed that such functional is not stable; performing high 
precision Hartree-Fock calculations (HF) using a sufficient 
number of iterations, the code converges to a solution 
where protons and neutrons are separated in clusters. The 
same conclusion was found for other Skyrme functionals, 
for example the LNS [11 one. Such unphysical configura- 
tions were discovered to be related to an excessive contri- 
bution of the term proportional to pi Api in the functional. 

Performing extensive HF calculations for even-even [2 
and odd-even nuclei |13) , we have found that there exists a 
relation between the position of the poles of the response 
function in SNM and the instabilities in finite nuclei. To 
guarantee that a nucleus will not manifest finite-size in- 
stabilities we have imposed the empirical constraint 

Ppolc > l-4/9sat, (4) 

where Psat is the saturation density of SNM. Such relation 
has to be satisfied for all the spin/isospin channels. The 
only instability that we allow to exist in the density region 
p < 1.4/9sat is the spinodal instability, that manifests in 
the (0,0,0) channel at low densities and low transferred 
momenta and characterizes the liquid-gas phase transition 
of nuclear matter. Using such criterion allows to quickly 
test large sets of Skyrme functionals. This is a great 



advantage compared to self consistent cranking calcula- 
tions 12| and systematic quasi-particle blocking calcula- 
tions |14j used to reveal the instabilities in spin channels. 

Usually, the time-odd terms of the functional are not 
directly constrained to observables. Indirect constraints 
may come from the fact that some functionals have been 
constructed with the requirements that the time odd terms 
lead to Landau parameters satisfying the conditions 



for £ = 0, 1 and equivalent ones for and G^. These 
conditions are used to prevent the appearance of global 
spin and/or isospin instabilities, but do not prevent the 
formation of finite size instabilities due to the gradient 
terms such as s^Ast and (Vst)^ for isospin t = 0,1. 

An example of such spin instability has been found with 
the SLy5 [TB] functional. In Figure [5] we show the posi- 
tion of the poles of the response function in SNM. We 
notice that there is a pole in the channel (1,1,0) close 
to saturation density that does not respect the criterium 
given in Eq. Performing Time-Dependent-Hartree- 

Fock (TDHF) calculations, Fracasso et al. [TB] have shown 
that the functional SLy5 is indeed unstable against spin- 
dipole excitations. 

4. Results 

Using the results discussed in the previous section, we now 
present a new fitting protocol to build stable Skyrme func- 
tionals. We use a modern version of the Saclay-Lyon fit- 
ting protocol similar to the one presented in the recent ar- 
ticle of Washiyama et al. [T7] . On top of such protocol we 
add the constraint given by Eq ^ . This constraint is en- 
forced at each iteration of the fit and does not significantly 
change the execution time since, as already discussed in 
ref. [1] , since the expressions of the response functions are 
analytical. 

Using the SLy5 functional as a starting point, we de- 
cided to build a new functional with similar properties. 
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Fig. 2: (Colors online). Position of the poles in SNM for two different functional. The horizontal black dotted-dashed line 
represents the saturation density of the system, paat, while the red dotted line indicates lApsat and it represents the constraint 
given in EqU) See text for details. 



but without instabilities in the spin channel. The result- 
ing functional has been called SLy5* and its coupling con- 
stants are reported in Table [TJ In the same table we com- 
pare nuclear matter properties at the saturation point for 
the new SLy5* and the original SLy5, showing that they 
all are essentially equals but the Thomas- Reiche-Kuhn en- 
hancement factor Ky which is larger with SLy5*. The value 
obtained with SLy5* is in agreement with the empirical 
value EO] close to k^^^p « 0.4. 



Parameter 


SLy5* 


SLy5 


ta [MeVfm^] 


-2495.310 


-2484.880 


h [MeVfm^] 


484.020 


483.130 


t2 [MeVfm^] 


-469.480 


-549.400 


t3 [MeVfm3+3"] 


13867.430 


13763.000 


Xo 


0.620 


0.778 


X\ 


-0.086 


-0.328 




-0.947 


-1.000 




0.934 


1.267 


a 


1/6 


1/6 


W [MeVfrn^] 


120.250 


126.000 


E/A [MeV] 


-16.02 


-15.98 


Psat [fm^^] 


0.160 


0.160 


m* /m 


0.700 


0.697 


Koo [MeV] 


229.8 


229.9 


a I [MeV] 


32.01 


32.03 




0.418 


0.250 



Table 1 : Parameters of the Skyrme functionals used in the 
article. In the last lines we give the main SNM properties 
around saturation for these functionals: binding energy 
per nucleon E/A, saturation density psat, isoscalar effec- 
tive mass m* /m, incompressibility Kao, symmetry energy 
coefficient aj and the Thomas-Reiche-Kuhn enhancement 
factor Ky. 



4.1. Equations of State 

In order to compare the isoscalar and isovector properties 
of the functionals SLy5 and SLy5* we have calculated the 



equation of state (EoS) of SNM and PNM as well as the 
neutron effective mass in these to extreme media. 

In the upper panel of Figure [31 we show the equations 
of state for SNM and Pure Neutron Matter (PNM). We 
observe a very good agreement between the results given 
by two functionals. In SNM the two curves stay on top of 
each other and in PNM this is the case up to 2 to 3 times 
saturation density. 




Fig. 3: (Colors online) In the upper panels, we show the equa- 
tion of state in SNM and PNM for the two forces used in the 
article: SLy5 (left panel) and SLy5* (right panel). In the lower 
panels, the neutron effective mass is depicted for the two sys- 
tems. See text for details. 

In the lower panel of the same figure, we report the 
evolution of the neutron effective mass in SNM and PNM. 
In SNM, the curves that give the neutron effective mass 
m* /m as a function of the density stay on top of each and 
can not be distinguished. In the PNM case, the neutron 
effective mass obtained from SLy5 is systematically lower 
than the effective mass in SNM. 

One can also observe that SLy5* almost gives the same 
evolution for the neutron effective mass in SNM and in 
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PNM. The origin of this property can be clarified by writ- 
ing the expressions for the neutron effective masse in the 
two media in terms of couphng constants |5] 



2m* 



2m* 



2m 
2m 



CIp (SNM), (6) 

{CI^GDp (PNM). (7) 



These couphng constants are CJ — 55.18 MeVfm^ and 
CI = 1.23 MeVfm^ for the new SLy5* and = 
56.25 MeVfm^ and C[ = 23.95 MeVfm^ for the original 
SLy5. We observe that while the two isoscalar coupling 
constants are very similar the isovector one are very dif- 
ferent. This, according to Eq. ©, explains the difference 
between the effective masses predicted by SLy5 and SLy5* 
as seen in Figure |3l The fact that SLy5* gives similar ef- 
fective mass for the neutron in SNM and PNM is due to 
the very small value of obtained with this interaction. 

4.2. Binding energies 

To further test the quality of the new SLy5* interaction, 
we calculated the binding energies of a set of semi-magic 
nuclei using the Hartree-Fock-Bogoliubov (HFB) method. 
For the particle-hole channel we use the two functionals 
SLy5 and SLy5* while for the particle-particle channel we 
used a contact pairing interaction surface-type. 

With such choice we obtain an average pairing gap [15] 
Aavcr = 1-201 MeV for SLy5 and Aavor = 1-164 MeV for 
SLy5*. We compare the resulting binding energies with 
the experimental ones given by Audi et al. [TH]- We ob- 
serve that there is no significant qualitative difference be- 
tween the two functionals for this set of nuclei. The bind- 
ing energies predicted with SLy5* are in average slightly 
smaller than with SLy5 but both functional compare with 
the experimental data with the same average deviation. 

4.3. Charge radii 

In Figl5]we compare the charge radii obtained with the two 
functionals against the experimental results taken from 
the UNEDF web page0. Once again, we observe that SLy5 
and SLy5* show the same behavior and the agreement 
between theory and experiment is of the same level. 

4.4. Stability against spin polarization 

Following the work done in Ref . [12] , we tested the stability 
of the new force SLy5* in self-consistent cranking calcu- 
lations. For this purpose, we calculated a high-spin state 
of the ground super-deformed band of of ^^'*Hg, that is 
at Jz = 54 (natural units), going up to several thousands 
iterations. We observe that SLy5* is perfectly stable and 
does not display the problems detected with SLy5, which 
typically leads to systematic non-convergence within few 
iterations. 

5. Conclusions 
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Fig. 4: (Colors onIine)In the upper panel, we show the dif- 
ference in biding energies for different isotopic chains of semi- 
magic nuclei obtained with SLy5 and SLy5*. In the lower panel 
the same, but for isotonic chains. 
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We have shown that, with the inclusion of the LR method 
into the fitting protocol, it is possible to build functionals 
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Fig. 5: (Colors online) Comparing charge radii for two different 
force for four different isotopic chains. 
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that do not present instabilities. As a first example, we 
used this improved protocol to cure the spin instability 
present in SLy5. The resulting functional has a compara- 
ble level of accuracy as the original one for the prediction 
of relevant observables as binding energies and charge radii 
for the set of semi-magic nuclei considered here, moreover 
most of the relevant infinite nuclear matter properties are 
unchanged only the Thomas-Reiche-Kuhn enhancement 
factor is significantly larger. 

Such method clearly demonstrates that it is possible to 
use the LR method for additional constraints to build a 
functional with no finite size instabilities in infinite nuclear 
matter and nuclei. 
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